Birthday data case study: Diagnosis of the HSGP approximation

Load packages

library(rstan)
library(cmdstanr)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = FALSE)
library(posterior)
library(bayesplot)
library(ggplot2)
library(loo)
library(latex2exp)
# library(reshape2)
# library(plyr)

Data

data <- read.csv(file="births_usa_1969.csv")
str(data)

'data.frame':   7305 obs. of  8 variables:
 $ year        : int  1969 1969 1969 1969 1969 1969 1969 1969 1969 1969 ...
 $ month       : int  1 1 1 1 1 1 1 1 1 1 ...
 $ day         : int  1 2 3 4 5 6 7 8 9 10 ...
 $ births      : int  8486 9002 9542 8960 8390 9560 9738 9734 9434 10042 ...
 $ day_of_year : int  1 2 3 4 5 6 7 8 9 10 ...
 $ day_of_week : int  3 4 5 6 7 1 2 3 4 5 ...
 $ id          : int  1 2 3 4 5 6 7 8 9 10 ...
 $ day_of_year2: int  1 2 3 4 5 6 7 8 9 10 ...

Input vector

x <- scale(data$id[], center=TRUE, scale=TRUE)

Output vector

y <- scale(data$births[], center=TRUE, scale=TRUE)

Standard deviation and mean of the output

std_y <- attr(y,"scaled:scale")
m_y <- attr(y,"scaled:center")

Standard deviation of the input

std_x <- attr(x,"scaled:scale")
std_x

[1] 2108.916

Year and week periods

period_year <- 365.25/std_x
period_week <- 7/std_x

Day of the year2 is a variable that indicates the correspondence between all the days with the days of a leap-year

day_of_year2 <- data$day_of_year2

Half range of the input domain

S <- (max(x) - min(x))/2
S

[1] 1.731695

Relationships among \(m\), \(l\) and \(c\)

1. Relationship among the number of basis functions \(m\), the lengthscale \(l\) and the boundary factor \(c\) for a squared exponential kernel:

\[ m= 1.73 \, \frac{c}{\left(\frac{l}{S}\right)^{1.05}}, \;\;\; {\small\text{or equivalently,}} \;\;\; \frac{l}{S}= 1.68 \, \frac{c^{0.95}}{m^{0.95}}\]

with \(c \, \geq \, 0.97 + 1.45 \, \frac{l}{S}\), and \(S\) being the half range of the input domain.

2. Relationship between \(m\) and \(l\) for a squared exponential periodic kernel:

\[m= 3.31 \, \frac{1}{l^{1.05}}, \;\;\; {\small\text{or equivalently,}} \;\;\; l= 3.12 \, \frac{1}{m^{0.95}}\]

Model

Let \(y_t\) denote the number of births on the \(t\)’th day,

\[\begin{align*} &y_{t} \sim \text{Normal}(\mu(t),\sigma^2), \\ &\mu(t) = f_1(t) + f_2(t) + f_3(t) + f_4(t) \end{align*}\]

The component \(f_1(t)\) represents the long-term trends modeled by a GP with squared exponential covariance function,

\[\begin{equation*} f_1(t) \sim \text{GP}(0,k_1), \hspace{5mm} k_1(t,t') = \alpha_1 \exp\!\!\left(\!-\frac{1}{2} \frac{(t-t')^2}{\ell_1^2}\right) \end{equation*}\]

The component \(f_2(t)\) represents the yearly smooth seasonal pattern, using a periodic squared exponential covariance function (with period 365.25 to match the average length of the year) in a GP model,

\[\begin{align*} &f_2(t) \sim \text{GP}(0,k_2), \\ &k_2(t,t') = \alpha_2 \exp\!\!\left(\!-\frac{2\,\sin^{\!2}\!(\pi(t-t')/365.25}{\ell_2^2}\right) \end{align*}\]

The component \(f_3(t)\) represents the weekly smooth pattern using a periodic squared exponential covariance function (with period 7 of length of the week) in a GP model,

\[\begin{align*} &f_3(t) \sim \text{GP}(0,k_3), \\ &k_3(t,t') = \alpha_3 \exp\!\!\left(\!-\frac{2\,\sin^{\!2}\!(\pi(t-t')/7}{\ell_3^2}\right) \end{align*}\]

The component \(f_4(t)\) represents the special days effects modeled as t-student distribution:

\[\begin{align*} &f_4(t) \sim \text{t-student}(0,\alpha_4^2), \; {\small\text{with 1 degree of freedom}},\\ &\alpha_4 \sim \text{Normal}(0,0.1) \end{align*}\]

Model in Stan code

//saved in 'stancode_def.stan'
functions {
  vector diagSPD_EQ(real gpscale, real lscale, real L, int M) {
    return sqrt((gpscale^2) * sqrt(2*pi()) * lscale * exp(-0.5*(lscale*pi()/2/L)^2 * linspaced_vector(M, 1, M)^2));
  }
  vector diagSPD_periodic(real gpscale, real lscale, int M) {
    real a = 1/lscale^2;
    int one_to_M[M];
    for (m in 1:M) one_to_M[m] = m;
    vector[M] q = sqrt(gpscale^2 * 2 / exp(a) * to_vector(modified_bessel_first_kind(one_to_M, a)));
    return append_row(q,q);
  }
  matrix PHI_EQ(int N, int M, real L, vector x) {
     matrix[N,M] PHI = sin(diag_post_multiply(rep_matrix(pi()/(2*L) * (x+L), M), linspaced_vector(M, 1, M)))/sqrt(L);
     for (m in 1:M)
       PHI[,m] = PHI[,m] - mean(PHI[,m]);
     return PHI;
  }
  matrix PHI_periodic(int N, int M, real w0, vector x) {
    matrix[N,M] mw0x = diag_post_multiply(rep_matrix(w0*x, M), linspaced_vector(M, 1, M));
    matrix[N,M] PHI = append_col(cos(mw0x), sin(mw0x));
    for (m in 1:M)
      PHI[,m] = PHI[,m] - mean(PHI[,m]);
    return PHI;
  }
}
data {
    real c_f1;                  //boundary value for function 1
    int<lower=1> M_f1;          //num basis functions for function 1
    int<lower=-1> J_f2;         //num cosine and sinu functions for function 3
    int<lower=-1> J_f3;         //num cosine and sinu functions for function 4
    int<lower=1> N;             //num observations
    vector[N] x;                //input vector
    vector[N] y;                //target vector
    real period_year;           //period of the year
    real period_week;           //period of the week
    int day_of_year2[N];        //day of the year inside a leap-year
}
transformed data {
    real L_f1= c_f1*max(fabs(x));
    //Basis functions for f1, f2 and f3
    matrix[N,M_f1] PHI_f1 = PHI_EQ(N, M_f1, L_f1, x);
    matrix[N,2*J_f2] PHI_f2 = PHI_periodic(N, J_f2, 2*pi()/period_year, x);
    matrix[N,2*J_f3] PHI_f3 = PHI_periodic(N, J_f3, 2*pi()/period_week, x);
}
parameters {
    vector[N] f;
    real intercept;
    //variables for the basis function models
    vector[M_f1] beta_f1;
    vector[2*J_f2] beta_f2;
    vector[2*J_f3] beta_f3;
    //hyperparameters
    vector<lower=0>[3] lscale;
    vector<lower=0>[3] gpscale;
    real<lower=0> noise;
    //
    //t-student prior for special days effects
    vector[366] f4;
    real<lower=0> sigma_f4;
}
transformed parameters{
    vector[N] f1;
    vector[N] f2;
    vector[N] f3;
    {
    vector[M_f1] diagSPD_f1 = diagSPD_EQ(gpscale[1], lscale[1], L_f1, M_f1);
    vector[2*J_f2] diagSPD_f2 = diagSPD_periodic(gpscale[2], lscale[2], J_f2);
    vector[2*J_f3] diagSPD_f3 = diagSPD_periodic(gpscale[3], lscale[3], J_f3);
    //
    vector[M_f1] SPD_beta_f1 = diagSPD_f1 .* beta_f1;
    vector[2*J_f2] SPD_beta_f2 = diagSPD_f2 .* beta_f2;
    vector[2*J_f3] SPD_beta_f3 = diagSPD_f3 .* beta_f3;
    //
    f1 = PHI_f1[,] * SPD_beta_f1;
    f2 = PHI_f2[,] * SPD_beta_f2;
    f3 = PHI_f3[,] * SPD_beta_f3;
    }
}
model{
    intercept ~ normal(0,1);
    beta_f1 ~ normal(0,1);
    beta_f2 ~ normal(0,1);
    beta_f3 ~ normal(0,1);
    //
    lscale ~ normal(0,2);           //GP lengthscales
    gpscale ~ normal(0,10);         //GP magnitudes
    noise ~ normal(0,1);            //model noise
    //
    //t-student prior for special days effects
    f4 ~ student_t(1, 0, sigma_f4);
    sigma_f4 ~ normal(0, 0.1);
    //
    f= intercept + f1 + f2 + f3 + f4[day_of_year2];
    target += normal_lpdf(y | f, noise);
}
generated quantities{
  vector[N] y_rep;
  vector[N] log_lik;
  for(n in 1:N){
    y_rep[n] = normal_rng(f[n], noise);
    log_lik[n] = normal_lpdf(y[n] | f[n], noise);
  }
}

User-guide for diagnosis

Assumption of the diagnosis tool:

• Under inaccurate HSGP approximation, the estimated lengthscale \(\hat{l}\) is (always) smaller than the true lengthscale \(l\).

User-guide with the steps to perform diagnosis:

0. Make a first guess of the lengthscale \(l^1\) of the function to be learned.

Iteration 1

1. Obtain the minimum valid boundary factor \(c^1\) determined by the first guess \(l^1\), by using the functional relationship \(\, c^1 \, \geq \, 0.97 + 1.45 \, \frac{l^1}{S}\), where \(S\) represents the half range of the input domain.

2. Obtain the mimimum valid number of basis functions \(m^1\) determined by \(l^1\) and \(c^1\), by using the functional relationship \(\, m^1= 1.73 \, \frac{c^1}{\left(\frac{l^1}{S}\right)^{1.05}}\). Notice that \(l^1\) can also be read as the minimum lengthscale that can be accurately fitted determined by \(m^1\) and \(c^1\) \(\left(\frac{l^1}{S}= 1.68 \, \frac{(c^1)^{0.95}}{(m^1)^{0.95}}\right)\).

3. Fit the HSGP model and assess residuals: check any residual trends, compute rmse (Root mean square error), \(R^2\) (Coefficient of variation) and elpd (Expected log predicitve density).

4. Check whether the estimated \(\hat{l^1}\) is equal to or greater than that \(l^1\) (the minimum lengthscale that can be accurately fitted determined by the \(m^1\) and \(c^1\) used).

5. If the verification in step 4 is TRUE, the HSGP model approximation must be sufficiently accurate, and diagnosis ends here. Otherwise, continue with the next step 6.

Iteration 2

6. Use the estimated lengthscale \(\hat{l^1}\) to obtain the new \(c^2\).

7. Repeat steps 1-5 and update parameters.

Note: In order to be a bit more conservative and make as less iterations as possible to get an accurate approximation, 5 or 10 basis functions can be added to the computed \(m\) in each iteration.

Fitting and diagnosis

Initializing programming objects:

standata <- list()
fit <- list()
l_f1 <- vector()
c_f1 <- vector()
m_f1 <- vector()
l_f2 <- vector()
m_f2 <- vector()
l_f3 <- vector()
m_f3 <- vector()
l_f1_hat <- vector()
l_f2_hat <- vector()
l_f3_hat <- vector()
res <- list()
diagnosis <- list()
rmse <- vector()
eR2 <- vector()
elpd <- vector()

Iteration 1

Iteration index

i <- 1

Setting \(m\), \(l\) and \(c\)

Smooth trend function \(f_1\)

1. Making the first guess that the lengthscale \(l\) might correspond to around 3 years (3*365=1095 days) of the input dimension

l_f1[i] <- round(1095 / std_x, 2)
l_f1[i]

[1] 0.52

2. The boundary factor \(c\) has to fulfill \(c \, \geq \, 0.97 + 1.45 \, \frac{l}{S}\)

0.97 + 1.45*l_f1[i]/S

[1] 1.405412

Then, \(c\) is set to

c_f1[i] <- 1.5
c_f1[i]

[1] 1.5

3. The number of basis functions \(m\) as a function of \(l\): \(\; m= 1.73 \, \frac{c}{(l/S)^{1.05}}\)

m_f1[i] <- ceiling(1.73 * c_f1[i]/(l_f1[i]/S)^1.05)
m_f1[i]

[1] 10

Yearly periodic effects function \(f_2\)

1. Making the first guess that the lengthscale \(l\) might corresponds to half of the period of \(f_2\). So \(l=0.5\)

l_f2[i] <- 0.5
l_f2[i]

[1] 0.5

2. Number of basis functions \(m\) as a function of \(l\): \(\; m= 3.31 \, \frac{1}{l^{1.05}}\)

m_f2[i] <-  ceiling(3.31 * 1/l_f2[i]^1.05)
m_f2[i]

[1] 7

Weekly periodic effects function \(f_3\)

1. Making the first guess that the lengthscale \(l\) might corresponds to half of the period of \(f_3\). So \(l=0.5\)

l_f3[i] <- 0.5
l_f3[i]

[1] 0.5

2. Number of basis functions \(m\) as a function of \(l\): \(\; m= 3.31 \, \frac{1}{l^{1.05}}\)

m_f3[i] <- ceiling(3.31 * 1/l_f3[i]^1.05)
m_f3[i]

[1] 7

Data to Stan

standata[[i]] <- list(M_f1= m_f1[i], 
                         c_f1= c_f1[i], 
                         J_f2= m_f2[i], 
                         J_f3= m_f3[i], 
                         x= x[,1], 
                         y= y[,1], 
                         N= length(x), 
                         period_year= period_year, 
                         period_week= period_week,
                         day_of_year2= day_of_year2
)
str(standata[[i]])

List of 10
 $ M_f1        : num 10
 $ c_f1        : num 1.5
 $ J_f2        : num 7
 $ J_f3        : num 7
 $ x           : num [1:7305] -1.73 -1.73 -1.73 -1.73 -1.73 ...
 $ y           : num [1:7305] -1.0316 -0.5739 -0.0949 -0.6111 -1.1168 ...
 $ N           : int 7305
 $ period_year : num 0.173
 $ period_week : num 0.00332
 $ day_of_year2: int [1:7305] 1 2 3 4 5 6 7 8 9 10 ...

Model fitting

Compiling the model

# birthday_mod <- cmdstanr::cmdstan_model(stan_file = "stancode_def.stan")

Model sampling

load("fit.rData")

# fit[[i]] <- birthday_mod$sample(data= standata[[i]], iter_warmup=200, iter_sampling=200, chains=4, thin=4, init=0.5, adapt_delta=0.9, save_warmup=FALSE)
# fit[[i]] <- read_stan_csv(fit[[i]]$output_files())

# save(fit, file="fit.rData")

Summaries of variable estimates

param = c("intercept","lscale[1]","lscale[2]","lscale[3]","gpscale[1]","gpscale[2]","gpscale[3]","noise","sigma_f4")
summary(fit[[i]], pars = param, probs = c(0.025, 0.5, 0.975), digits_summary = 4)$summary

                 mean      se_mean          sd        2.5%        50%
intercept  0.01621316 0.0004436968 0.005327212 0.005182119 0.01648965
lscale[1]  0.24056074 0.0121119972 0.153241311 0.008576415 0.23887050
lscale[2]  0.43967898 0.0105689537 0.082971748 0.221308300 0.44917800
lscale[3]  1.58862805 0.0655106390 0.613293463 0.608230400 1.48990500
gpscale[1] 2.80090466 0.1469704409 1.946717801 1.290186500 2.18231500
gpscale[2] 0.39863947 0.0113248627 0.117024120 0.236473275 0.38862950
gpscale[3] 2.81596814 0.1902765926 2.358660650 0.745981300 1.97060500
noise      0.30919913 0.0001872498 0.002675033 0.304286125 0.30911650
sigma_f4   0.02119854 0.0003020931 0.003114659 0.016262200 0.02068870
                97.5%     n_eff      Rhat
intercept  0.02563535 144.15422 0.9931835
lscale[1]  0.49549075 160.07378 0.9885406
lscale[2]  0.57091045  61.63061 1.0446003
lscale[3]  2.94770575  87.64215 1.0328162
gpscale[1] 7.51044400 175.44702 1.0021600
gpscale[2] 0.68352850 106.77876 1.0550956
gpscale[3] 8.75818300 153.65977 1.0039020
noise      0.31408470 204.08715 1.0019957
sigma_f4   0.02781466 106.30151 1.0183450

Simulation chains for the variables after warmup

traceplot(fit[[i]], pars = param, include = TRUE, unconstrain = FALSE, inc_warmup = FALSE, window = NULL, nrow = NULL, ncol = NULL)

Plot of the mean posterior function \(f_1\)

Plot of the mean posterior function \(f_2\)

Plot of the mean posterior function \(f_1 + f_2\)

Plot of the residuals

Model evaluation

Residuals

f <- summary(fit[[i]], pars = c("f"))$summary[,1]
res[[i]] <- standata[[i]]$y - f

ggplot(as.data.frame(res[[i]]), aes(res[[i]]))  +
  geom_histogram(color = 'white') +
  theme_classic()

Root mean squared error

rmse[i] <- sqrt(mean(res[[i]]^2))
rmse[i]

[1] 0.3055324

Bayesian \(R^2\) (Coefficient of determination)

noise <- as.matrix(fit[[i]], pars = c("noise"))
sd_f <- apply(as.matrix(fit[[i]], pars = c("f")), 1, sd)
R2 <- sd_f^2/(sd_f^2 + noise^2)

ggplot(as.data.frame(R2), aes(R2))  +
  geom_histogram(color = 'white') +
  theme_classic()

Median of \(R^2\)

eR2[i] <- median(R2)
eR2[i]

[1] 0.904461

Log predictive density (lpd)

lpd <- summary(fit[[i]], pars = c("log_lik"))$summary[,1]

ggplot(as.data.frame(lpd), aes(lpd))  +
  geom_histogram(color = 'white') +
  theme_classic()

Median of lpd

elpd[i] <- median(lpd)
elpd[i]

[1] 0.07947658

Diagnosis

Smooth trend function \(f_1\)

1. Estimated lengthscale \(\hat{l}\)

l_f1_hat[i] <- round(summary(fit[[i]], pars = "lscale[1]")$summary[,1], 2)
l_f1_hat[i]

[1] 0.24

2. Check whether \(\hat{l}\) is equal to or greater than l_f1 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f1_hat[i] >= l_f1[i]

[1] FALSE

Yearly periodic effects function \(f_2\)

1. Estimated lengthscale \(\hat{l}\)

l_f2_hat[i] <- round(summary(fit[[i]], pars = "lscale[2]")$summary[,1], 2)
l_f2_hat[i]

[1] 0.44

2. Check whether \(\hat{l}\) is equal to or greater than l_f2 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f2_hat[i] >= l_f2[i]

[1] FALSE

Weekly periodic effects function \(f_3\)

1. Estimated lengthscale \(\hat{l}\)

l_f3_hat[i] <- round(summary(fit[[i]], pars = "lscale[3]")$summary[,1], 2)
l_f3_hat[i]

[1] 1.59

2. Check whether \(\hat{l}\) is equal to or greater than l_f3 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f3_hat[i] >= l_f3[i]

[1] TRUE

Summary table

iter	GP_func	l	c	m	l_hat	l_hat > l	rmse	R2	elpd
1	f1	0.52	1.5	10	0.24	FALSE	0.306	0.904	0.079
1	f2	0.50	NA	7	0.44	FALSE	0.306	0.904	0.079
1	f3	0.50	NA	7	1.59	TRUE	0.306	0.904	0.079

Iteration 2

Iteration index

i <- 2

Setting \(m\), \(l\) and \(c\)

Smooth trend function \(f_1\)

1. Updating \(l\) with its estimate \(\hat{l}\) in the first iteration

l_f1[i] <- l_f1_hat[i-1]
l_f1[i]

[1] 0.24

2. The boundary factor \(c\) has to fulfill \(c \, \geq \, 0.97 + 1.45 \, \frac{l}{S}\)

0.97 + 1.45*l_f1[i]/S

[1] 1.170959

Then, \(c\) is set to

c_f1[i] <- 1.5
c_f1[i]

[1] 1.5

3. The number of basis functions \(m\) as a function of \(l\): \(\; m= 1.73 \, \frac{c}{(l/S)^{1.05}}\)

if(1.73 * c_f1[i]/(l_f1[i]/S)^1.05 < 5) m_f1[i] <- 5 else m_f1[i] <- ceiling(1.73 * c_f1[i]/(l_f1[i]/S)^1.05)
m_f1[i]

[1] 21

Yearly periodic effects function \(f_2\)

1. Updating \(l\) with its estimate \(\hat{l}\) in the first iteration

l_f2[i] <- l_f2_hat[i-1]
l_f2[i]

[1] 0.44

2. Number of basis functions \(m\) as a function of \(l\): \(\; m= 3.31 \, \frac{1}{l^{1.05}}\)

if(3.31 * 1/l_f2[i]^1.05 < 5) m_f2[i] <- 5 else m_f2[i] <- ceiling( 3.31 * 1/l_f2[i]^1.05)
m_f2[i]

[1] 8

Weekly periodic effects function \(f_3\)

1. Updating \(l\) with its estimate \(\hat{l}\) in the first iteration

l_f3[i] <- l_f3_hat[i-1]
l_f3[i]

[1] 1.59

2. Number of basis functions \(m\) as a function of \(l\): \(\; m= 3.31 \, \frac{1}{l^{1.05}}\)

if(3.31 * 1/l_f3[i]^1.05 < 5) m_f3[i] <- 5 else m_f3[i] <- ceiling(3.31 * 1/l_f3[i]^1.05)
m_f3[i]

[1] 5

Data to Stan

standata[[i]] <- list(M_f1= m_f1[i], 
                 c_f1= c_f1[i], 
                 J_f2= m_f2[i], 
                 J_f3= m_f3[i], 
                 x= x[,1], 
                 y= y[,1], 
                 N= length(x), 
                 period_year= period_year, 
                 period_week= period_week,
                 day_of_year2= day_of_year2
)
str(standata[[i]])

List of 10
 $ M_f1        : num 21
 $ c_f1        : num 1.5
 $ J_f2        : num 8
 $ J_f3        : num 5
 $ x           : num [1:7305] -1.73 -1.73 -1.73 -1.73 -1.73 ...
 $ y           : num [1:7305] -1.0316 -0.5739 -0.0949 -0.6111 -1.1168 ...
 $ N           : int 7305
 $ period_year : num 0.173
 $ period_week : num 0.00332
 $ day_of_year2: int [1:7305] 1 2 3 4 5 6 7 8 9 10 ...

Model fitting

Compiling the model

# birthday_mod <- cmdstanr::cmdstan_model(stan_file = "stancode_def.stan")

Model sampling (using cmdstanr package)

# load("fit.rData")

# fit[[i]] <- birthday_mod$sample(data= standata[[i]], iter_warmup=200, iter_sampling=200, chains=4, thin=4, init=0.5, adapt_delta=0.9, save_warmup=FALSE)
# fit[[i]] <- rstan::read_stan_csv(fit[[i]]$output_files())

# save(fit, file="fit.rData")

Summaries of variable estimates

param = c("intercept","lscale[1]","lscale[2]","lscale[3]","gpscale[1]","gpscale[2]","gpscale[3]","noise","sigma_f4")
summary(fit[[i]], pars = param, probs = c(0.025, 0.5, 0.975), digits_summary = 4)$summary

                 mean      se_mean          sd       2.5%       50%       97.5%
intercept  0.03209490 0.0006963903 0.005952288 0.02002443 0.0324404  0.04281516
lscale[1]  0.20129869 0.0043869778 0.045159338 0.07162255 0.2089775  0.26623510
lscale[2]  0.34841957 0.0055318075 0.054283476 0.21823128 0.3536220  0.44857008
lscale[3]  1.70160588 0.0482643667 0.641194415 0.70088583 1.5790850  3.19423950
gpscale[1] 0.84056559 0.0285485120 0.330861998 0.44415168 0.7609500  1.68403925
gpscale[2] 0.34160833 0.0096154126 0.092702693 0.21751665 0.3231235  0.60874603
gpscale[3] 3.14635309 0.1900064632 2.570161204 0.77181205 2.2757550 10.69412250
noise      0.29952692 0.0001630606 0.002453281 0.29424758 0.2996545  0.30394945
sigma_f4   0.02045204 0.0004095576 0.002985429 0.01490792 0.0203633  0.02544959
               n_eff      Rhat
intercept   73.05710 1.0379945
lscale[1]  105.96545 0.9947814
lscale[2]   96.29443 1.0078471
lscale[3]  176.49249 1.0047101
gpscale[1] 134.31571 1.0031222
gpscale[2]  92.94988 1.0369656
gpscale[3] 182.97172 1.0002500
noise      226.35816 1.0050146
sigma_f4    53.13536 1.0398422

Simulation chains for the variables after warmup

traceplot(fit[[i]], pars = param, include = TRUE, unconstrain = FALSE, inc_warmup = FALSE, window = NULL, nrow = NULL, ncol = NULL)

Plot of the mean posterior function \(f_1\)

Plot of the mean posterior function \(f_2\)

Plot of the mean posterior function \(f_1 + f_2\)

Plot of the residuals

Model evaluation

Residuals

f <- summary(fit[[i]], pars = c("f"))$summary[,1]
res[[i]] <- standata[[i]]$y - f

ggplot(as.data.frame(res[[i]]), aes(res[[i]]))  +
  geom_histogram(color = 'white') +
  theme_classic()

Root mean squared error

rmse[i] <- sqrt(mean(res[[i]]^2))
rmse[i]

[1] 0.2954186

Bayesian \(R^2\) (Coefficient of determination)

noise <- as.matrix(fit[[i]], pars = c("noise"))
sd_f <- apply(as.matrix(fit[[i]], pars = c("f")), 1, sd)
R2 <- sd_f^2/(sd_f^2 + noise^2)

ggplot(as.data.frame(R2), aes(R2))  +
  geom_histogram(color = 'white') +
  theme_classic()

Median of \(R^2\)

eR2[i] <- median(R2)
eR2[i]

[1] 0.9101999

Log predictive density (lpd)

lpd <- summary(fit[[i]], pars = c("log_lik"))$summary[,1]

ggplot(as.data.frame(lpd), aes(lpd))  +
  geom_histogram(color = 'white') +
  theme_classic()

Median of lpd

elpd[i] <- median(lpd)
elpd[i]

[1] 0.1183452

Diagnosis

Smooth trend function \(f_1\)

1. Estimated lengthscale \(\hat{l}\)

l_f1_hat[i] <- round(summary(fit[[i]], pars = "lscale[1]")$summary[,1], 2)
l_f1_hat[i]

[1] 0.2

2. Check whether \(\hat{l}\) is equal to or greater than l_f1 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f1_hat[i] >= l_f1[i]

[1] FALSE

Yearly periodic effects function \(f_2\)

1. Estimated lengthscale \(\hat{l}\)

l_f2_hat[i] <- round(summary(fit[[i]], pars = "lscale[2]")$summary[,1], 2)
l_f2_hat[i]

[1] 0.35

2. Check whether \(\hat{l}\) is equal to or greater than l_f2 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f2_hat[i] >= l_f2[i]

[1] FALSE

Weekly periodic effects function \(f_3\)

1. Estimated lengthscale \(\hat{l}\)

l_f3_hat[i] <- round(summary(fit[[i]], pars = "lscale[3]")$summary[,1], 2)
l_f3_hat[i]

[1] 1.7

2. Check whether \(\hat{l}\) is equal to or greater than l_f3 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f3_hat[i] >= l_f3[i]

[1] TRUE

Summary table

iter	GP_func	l	c	m	l_hat	l_hat > l	rmse	R2	elpd
2	f1	0.24	1.5	21	0.20	FALSE	0.295	0.91	0.118
2	f2	0.44	NA	8	0.35	FALSE	0.295	0.91	0.118
2	f3	1.59	NA	5	1.70	TRUE	0.295	0.91	0.118

Iteration 3

Iteration index

i <- 3

Setting \(m\), \(l\) and \(c\)

Smooth trend function \(f_1\)

1. Updating \(l\) with its estimate \(\hat{l}\) in the second iteration

l_f1[i] <- l_f1_hat[i-1]
l_f1[i]

[1] 0.2

2. The boundary factor \(c\) has to fulfill \(c \, \geq \, 0.97 + 1.45 \, \frac{l}{S}\)

0.97 + 1.45*l_f1[i]/S

[1] 1.137466

Then, \(c\) is set to

c_f1[i] <- 1.5
c_f1[i]

[1] 1.5

3. The number of basis functions \(m\) as a function of \(l\): \(\; m= 1.73 \, \frac{c}{(l/S)^{1.05}}\)

if(1.73 * c_f1[i]/(l_f1[i]/S)^1.05 < 5) m_f1[i] <- 5 else m_f1[i] <- ceiling(1.73 * c_f1[i]/(l_f1[i]/S)^1.05)
m_f1[i]

[1] 26

Yearly periodic effects function \(f_2\)

1. Updating \(l\) with its estimate \(\hat{l}\) in the second iteration

l_f2[i] <- l_f2_hat[i-1]
l_f2[i]

[1] 0.35

2. Number of basis functions \(m\) as a function of \(l\): \(\; m= 3.31 \, \frac{1}{l^{1.05}}\)

if(3.31 * 1/l_f2[i]^1.05 < 5) m_f2[i] <- 5 else m_f2[i] <- ceiling(3.31 * 1/l_f2[i]^1.05)
m_f2[i]

[1] 10

Weekly periodic effects function \(f_3\)

1. Updating \(l\) with its estimate \(\hat{l}\) in the second iteration

l_f3[i] <- l_f3_hat[i-1]
l_f3[i]

[1] 1.7

2. Number of basis functions \(m\) as a function of \(l\): \(\; m= 3.31 \, \frac{1}{l^{1.05}}\)

if(3.31 * 1/l_f3[i]^1.05 < 5) m_f3[i] <- 5 else m_f3[i] <- ceiling(3.31 * 1/l_f3[i]^1.05)
m_f3[i]

[1] 5

Data to Stan

standata[[i]] <- list(M_f1= m_f1[i], 
                 c_f1= c_f1[i], 
                 J_f2= m_f2[i], 
                 J_f3= m_f3[i], 
                 x= x[,1], 
                 y= y[,1], 
                 N= length(x), 
                 period_year= period_year, 
                 period_week= period_week,
                 day_of_year2= day_of_year2
)
str(standata[[i]])

List of 10
 $ M_f1        : num 26
 $ c_f1        : num 1.5
 $ J_f2        : num 10
 $ J_f3        : num 5
 $ x           : num [1:7305] -1.73 -1.73 -1.73 -1.73 -1.73 ...
 $ y           : num [1:7305] -1.0316 -0.5739 -0.0949 -0.6111 -1.1168 ...
 $ N           : int 7305
 $ period_year : num 0.173
 $ period_week : num 0.00332
 $ day_of_year2: int [1:7305] 1 2 3 4 5 6 7 8 9 10 ...

Model fitting

Compiling the model

# birthday_mod <- cmdstanr::cmdstan_model(stan_file = "stancode_def.stan")

Model sampling

# load("fit.rData")

# fit[[i]] <- birthday_mod$sample(data= standata[[i]], iter_warmup=200, iter_sampling=200, chains=4, thin=4, init=0.5, adapt_delta=0.9, save_warmup=FALSE)
# fit[[i]] <- rstan::read_stan_csv(fit[[i]]$output_files())

# save(fit, file="fit.rData")

Summaries of variable estimates

param = c("intercept","lscale[1]","lscale[2]","lscale[3]","gpscale[1]","gpscale[2]","gpscale[3]","noise","sigma_f4")
summary(fit[[i]], pars = param, probs = c(0.025, 0.5, 0.975), digits_summary = 4)$summary

                 mean      se_mean          sd       2.5%        50%      97.5%
intercept  0.02351247 0.0005943055 0.004968132 0.01406945 0.02367795 0.03283830
lscale[1]  0.13123223 0.0066514750 0.052912420 0.00864844 0.14233000 0.20619690
lscale[2]  0.29457215 0.0027172111 0.031269295 0.23831755 0.29434700 0.35713853
lscale[3]  1.57049730 0.0544457001 0.642421873 0.49599778 1.51369500 2.83534325
gpscale[1] 0.90058899 0.0908635474 0.566716049 0.53282395 0.75180300 2.20484175
gpscale[2] 0.32634141 0.0067091670 0.073006389 0.21944840 0.31186150 0.51495445
gpscale[3] 2.96833263 0.1920584568 2.474883866 0.69133465 2.05472500 8.99802925
noise      0.29703793 0.0001877659 0.002682057 0.29178015 0.29692300 0.30235395
sigma_f4   0.01585057 0.0003702678 0.002837476 0.01091628 0.01561165 0.02172783
               n_eff      Rhat
intercept   69.88221 1.0264599
lscale[1]   63.28187 1.0210061
lscale[2]  132.43093 1.0198239
lscale[3]  139.22379 1.0266769
gpscale[1]  38.90018 1.0952120
gpscale[2] 118.40895 1.0111913
gpscale[3] 166.05149 1.0048522
noise      204.03411 0.9948262
sigma_f4    58.72630 1.0220584

Simulation chains for the variables after warmup

traceplot(fit[[i]], pars = param, include = TRUE, unconstrain = FALSE, inc_warmup = FALSE, window = NULL, nrow = NULL, ncol = NULL)

Plot of the mean posterior function \(f_1\)

Plot of the mean posterior function \(f_2\)

Plot of the mean posterior function \(f_1 + f_2\)

Plot of the residuals

Model evaluation

Residuals

f <- summary(fit[[i]], pars = c("f"))$summary[,1]
res[[i]] <- standata[[i]]$y - f

ggplot(as.data.frame(res[[i]]), aes(res[[i]]))  +
  geom_histogram(color = 'white') +
  theme_classic()

Root mean squared error

rmse[i] <- sqrt(mean(res[[i]]^2))
rmse[i]

[1] 0.2932521

Bayesian \(R^2\) (Coefficient of determination)

noise <- as.matrix(fit[[i]], pars = c("noise"))
sd_f <- apply(as.matrix(fit[[i]], pars = c("f")), 1, sd)
R2 <- sd_f^2/(sd_f^2 + noise^2)

ggplot(as.data.frame(R2), aes(R2))  +
  geom_histogram(color = 'white') +
  theme_classic()

Median of \(R^2\)

eR2[i] <- median(R2)
eR2[i]

[1] 0.9118072

Log predictive density (lpd)

lpd <- summary(fit[[i]], pars = c("log_lik"))$summary[,1]

ggplot(as.data.frame(lpd), aes(lpd))  +
  geom_histogram(color = 'white') +
  theme_classic()

Median of lpd

elpd[i] <- median(lpd)
elpd[i]

[1] 0.1315323

Diagnosis

Smooth trend function \(f_1\)

1. Estimated lengthscale \(\hat{l}\)

l_f1_hat[i] <- round(summary(fit[[i]], pars = "lscale[1]")$summary[,1], 2)
l_f1_hat[i]

[1] 0.13

2. Check whether \(\hat{l}\) is equal to or greater than l_f1 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f1_hat[i] >= l_f1[i]

[1] FALSE

Yearly periodic effects function \(f_2\)

1. Estimated lengthscale \(\hat{l}\)

l_f2_hat[i] <- round(summary(fit[[i]], pars = "lscale[2]")$summary[,1], 2)
l_f2_hat[i]

[1] 0.29

2. Check whether \(\hat{l}\) is equal to or greater than l_f2 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f2_hat[i] >= l_f2[i]

[1] FALSE

Weekly periodic effects function \(f_3\)

1. Estimated lengthscale \(\hat{l}\)

l_f3_hat[i] <- round(summary(fit[[i]], pars = "lscale[3]")$summary[,1], 2)
l_f3_hat[i]

[1] 1.57

2. Check whether \(\hat{l}\) is equal to or greater than l_f3 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f3_hat[i] >= l_f3[i]

[1] FALSE

Summary table

iter	GP_func	l	c	m	l_hat	l_hat > l	rmse	R2	elpd
3	f1	0.20	1.5	26	0.13	FALSE	0.293	0.912	0.132
3	f2	0.35	NA	10	0.29	FALSE	0.293	0.912	0.132
3	f3	1.70	NA	5	1.57	FALSE	0.293	0.912	0.132

Iteration 4

Iteration index

i <- 4

Setting \(m\), \(l\) and \(c\)

Smooth trend function \(f_1\)

1. Updating \(l\) with its estimate \(\hat{l}\) in the third iteration

l_f1[i] <- l_f1_hat[i-1]
l_f1[i]

[1] 0.13

2. The boundary factor \(c\) has to fulfill \(c \, \geq \, 0.97 + 1.45 \, \frac{l}{S}\)

0.97 + 1.45*l_f1[i]/S

[1] 1.078853

Then, \(c\) is set to

c_f1[i] <- 1.5
c_f1[i]

[1] 1.5

3. The number of basis functions \(m\) as a function of \(l\): \(\; m= 1.73 \, \frac{c}{(l/S)^{1.05}}\)

if(1.73 * c_f1[i]/(l_f1[i]/S)^1.05 < 5) m_f1[i] <- 5 else m_f1[i] <- ceiling(1.73 * c_f1[i]/(l_f1[i]/S)^1.05)
m_f1[i]

[1] 40

Yearly periodic effects function \(f_2\)

1. Updating \(l\) with its estimate \(\hat{l}\) in the third iteration

l_f2[i] <- l_f2_hat[i-1]
l_f2[i]

[1] 0.29

2. Number of basis functions \(m\) as a function of \(l\): \(\; m= 3.31 \, \frac{1}{l^{1.05}}\)

if(3.31 * 1/l_f2[i]^1.05 < 5) m_f2[i] <- 5 else m_f2[i] <- ceiling(3.31 * 1/l_f2[i]^1.05)
m_f2[i]

[1] 13

Weekly periodic effects function \(f_3\)

1. Updating \(l\) with its estimate \(\hat{l}\) in the third iteration

l_f3[i] <- l_f3_hat[i-1]
l_f3[i]

[1] 1.57

2. Number of basis functions \(m\) as a function of \(l\): \(\; m= 3.31 \, \frac{1}{l^{1.05}}\)

if(3.31 * 1/l_f3[i]^1.05 < 5) m_f3[i] <- 5 else m_f3[i] <- ceiling(3.31 * 1/l_f3[i]^1.05)
m_f3[i]

[1] 5

Data to Stan

standata[[i]] <- list(M_f1= m_f1[i], 
                 c_f1= c_f1[i], 
                 J_f2= m_f2[i], 
                 J_f3= m_f3[i], 
                 x= x[,1], 
                 y= y[,1], 
                 N= length(x), 
                 period_year= period_year, 
                 period_week= period_week,
                 day_of_year2= day_of_year2
)
str(standata[[i]])

List of 10
 $ M_f1        : num 40
 $ c_f1        : num 1.5
 $ J_f2        : num 13
 $ J_f3        : num 5
 $ x           : num [1:7305] -1.73 -1.73 -1.73 -1.73 -1.73 ...
 $ y           : num [1:7305] -1.0316 -0.5739 -0.0949 -0.6111 -1.1168 ...
 $ N           : int 7305
 $ period_year : num 0.173
 $ period_week : num 0.00332
 $ day_of_year2: int [1:7305] 1 2 3 4 5 6 7 8 9 10 ...

Model fitting

Compiling the model

# birthday_mod <- cmdstanr::cmdstan_model(stan_file = "stancode_def.stan")

Model sampling

# load("fit.rData")

# fit[[i]] <- birthday_mod$sample(data= standata[[i]], iter_warmup=200, iter_sampling=200, chains=4, thin=4, init=0.5, adapt_delta=0.9, save_warmup=FALSE)
# fit[[i]] <- rstan::read_stan_csv(fit[[i]]$output_files())

# save(fit, file="fit.rData")

Summaries of variable estimates

param = c("intercept","lscale[1]","lscale[2]","lscale[3]","gpscale[1]","gpscale[2]","gpscale[3]","noise","sigma_f4")
summary(fit[[i]], pars = param, probs = c(0.025, 0.5, 0.975), digits_summary = 4)$summary

                 mean      se_mean          sd        2.5%       50%      97.5%
intercept  0.02163107 0.0003909226 0.004592181 0.011938308 0.0215220 0.03002273
lscale[1]  0.16195171 0.0010754892 0.013771676 0.133884550 0.1614695 0.19024127
lscale[2]  0.24561287 0.0024693364 0.028699315 0.181404650 0.2468775 0.29111728
lscale[3]  1.64165354 0.0552050194 0.636731523 0.516088400 1.5913050 2.98496650
gpscale[1] 0.62111567 0.0121793670 0.154971086 0.416687025 0.5857125 0.97306960
gpscale[2] 0.31635318 0.0067557442 0.061739874 0.218932425 0.3028695 0.47106488
gpscale[3] 3.01882210 0.2310811166 2.396089857 0.745482025 2.3385300 8.75506500
noise      0.29271492 0.0001838906 0.002583868 0.287375825 0.2928455 0.29797135
sigma_f4   0.01465976 0.0007576488 0.002767662 0.009466402 0.0144705 0.02018042
               n_eff      Rhat
intercept  137.99281 0.9992550
lscale[1]  163.96889 0.9869280
lscale[2]  135.07735 1.0386769
lscale[3]  133.03183 1.0004049
gpscale[1] 161.90189 0.9946525
gpscale[2]  83.51896 1.0225895
gpscale[3] 107.51704 1.0053584
noise      197.43379 1.0274909
sigma_f4    13.34413 1.2271513

Simulation chains for the variables after warmup

traceplot(fit[[i]], pars = param, include = TRUE, unconstrain = FALSE, inc_warmup = FALSE, window = NULL, nrow = NULL, ncol = NULL)

Plot of the mean posterior function \(f_1\)

Plot of the mean posterior function \(f_2\)

Plot of the mean posterior function \(f_1 + f_2\)

Plot of the residuals

Model evaluation

Residuals

f <- summary(fit[[i]], pars = c("f"))$summary[,1]
res[[i]] <- standata[[i]]$y - f

ggplot(as.data.frame(res[[i]]), aes(res[[i]]))  +
  geom_histogram(color = 'white') +
  theme_classic()

Root mean squared error

rmse[i] <- sqrt(mean(res[[i]]^2))
rmse[i]

[1] 0.2893248

Bayesian \(R^2\) (Coefficient of determination)

noise <- as.matrix(fit[[i]], pars = c("noise"))
sd_f <- apply(as.matrix(fit[[i]], pars = c("f")), 1, sd)
R2 <- sd_f^2/(sd_f^2 + noise^2)

ggplot(as.data.frame(R2), aes(R2))  +
  geom_histogram(color = 'white') +
  theme_classic()

Median of \(R^2\)

eR2[i] <- median(R2)
eR2[i]

[1] 0.9142275

Log predictive density (lpd)

lpd <- summary(fit[[i]], pars = c("log_lik"))$summary[,1]

ggplot(as.data.frame(lpd), aes(lpd))  +
  geom_histogram(color = 'white') +
  theme_classic()

Median of lpd

elpd[i] <- median(lpd)
elpd[i]

[1] 0.1487322

Diagnosis

Smooth trend function \(f_1\)

1. Estimated lengthscale \(\hat{l}\)

l_f1_hat[i] <- round(summary(fit[[i]], pars = "lscale[1]")$summary[,1], 2)
l_f1_hat[i]

[1] 0.16

2. Check whether \(\hat{l}\) is equal to or greater than l_f1 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f1_hat[i] >= l_f1[i]

[1] TRUE

Yearly periodic effects function \(f_2\)

1. Estimated lengthscale \(\hat{l}\)

l_f2_hat[i] <- round(summary(fit[[i]], pars = "lscale[2]")$summary[,1], 2)
l_f2_hat[i]

[1] 0.25

2. Check whether \(\hat{l}\) is equal to or greater than l_f2 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f2_hat[i] >= l_f2[i]

[1] FALSE

Weekly periodic effects function \(f_3\)

1. Estimated lengthscale \(\hat{l}\)

l_f3_hat[i] <- round(summary(fit[[i]], pars = "lscale[3]")$summary[,1], 2)
l_f3_hat[i]

[1] 1.64

2. Check whether \(\hat{l}\) is equal to or greater than l_f3 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f3_hat[i] >= l_f3[i]

[1] TRUE

Summary table

iter	GP_func	l	c	m	l_hat	l_hat > l	rmse	R2	elpd
4	f1	0.13	1.5	40	0.16	TRUE	0.289	0.914	0.149
4	f2	0.29	NA	13	0.25	FALSE	0.289	0.914	0.149
4	f3	1.57	NA	5	1.64	TRUE	0.289	0.914	0.149

Iteration 5

Iteration index

i <- 5

Setting \(m\), \(l\) and \(c\)

Smooth trend function \(f_1\)

1. Updating \(l\) with its estimate \(\hat{l}\) in the fourth iteration

l_f1[i] <- l_f1_hat[i-1]
l_f1[i]

[1] 0.16

2. The boundary factor \(c\) has to fulfill \(c \, \geq \, 0.97 + 1.45 \, \frac{l}{S}\)

0.97 + 1.45*l_f1[i]/S

[1] 1.103973

Then, \(c\) is set to

c_f1[i] <- 1.5
c_f1[i]

[1] 1.5

3. The number of basis functions \(m\) as a function of \(l\): \(\; m= 1.73 \, \frac{c}{(l/S)^{1.05}}\)

if(1.73 * c_f1[i]/(l_f1[i]/S)^1.05 < 5) m_f1[i] <- 5 else m_f1[i] <- ceiling(1.73 * c_f1[i]/(l_f1[i]/S)^1.05)

m_f1[i]

[1] 32

Yearly periodic effects function \(f_2\)

1. Updating \(l\) with its estimate \(\hat{l}\) in the third iteration

l_f2[i] <- l_f2_hat[i-1]
l_f2[i]

[1] 0.25

2. Number of basis functions \(m\) as a function of \(l\): \(\; m= 3.31 \, \frac{1}{l^{1.05}}\)

if(3.31 * 1/l_f2[i]^1.05 < 5) m_f2[i] <- 5 else m_f2[i] <- ceiling(3.31 * 1/l_f2[i]^1.05);
m_f2[i]

[1] 15

Weekly periodic effects function \(f_3\)

1. Updating \(l\) with its estimate \(\hat{l}\) in the fourth iteration

l_f3[i] <- l_f3_hat[i-1]
l_f3[i]

[1] 1.64

2. Number of basis functions \(m\) as a function of \(l\): \(\, m= 3.31 \, \frac{1}{l^{1.05}}\)

if(3.31 * 1/l_f3[i]^1.05 < 5) m_f3[i] <- 5 else m_f3[i] <- ceiling(3.31 * 1/l_f3[i]^1.05)
m_f3[i]

[1] 5

Data to Stan

standata[[i]] <- list(M_f1= m_f1[i], 
                 c_f1= c_f1[i], 
                 J_f2= m_f2[i], 
                 J_f3= m_f3[i], 
                 x= x[,1], 
                 y= y[,1], 
                 N= length(x), 
                 period_year= period_year, 
                 period_week= period_week,
                 day_of_year2= day_of_year2
)
str(standata[[i]])

List of 10
 $ M_f1        : num 32
 $ c_f1        : num 1.5
 $ J_f2        : num 15
 $ J_f3        : num 5
 $ x           : num [1:7305] -1.73 -1.73 -1.73 -1.73 -1.73 ...
 $ y           : num [1:7305] -1.0316 -0.5739 -0.0949 -0.6111 -1.1168 ...
 $ N           : int 7305
 $ period_year : num 0.173
 $ period_week : num 0.00332
 $ day_of_year2: int [1:7305] 1 2 3 4 5 6 7 8 9 10 ...

Model fitting

Compiling the model

# birthday_mod <- cmdstanr::cmdstan_model(stan_file = "stancode_def.stan")

Model sampling

# load("fit.rData")

# fit[[i]] <- birthday_mod$sample(data= standata[[i]], iter_warmup=200, iter_sampling=200, chains=4, thin=4, init=0.5, adapt_delta=0.9, save_warmup=FALSE)
# fit[[i]] <- rstan::read_stan_csv(fit[[i]]$output_files())

# save(fit, file="fit.rData")

Summaries of variable estimates

param = c("intercept","lscale[1]","lscale[2]","lscale[3]","gpscale[1]","gpscale[2]","gpscale[3]","noise","sigma_f4")
summary(fit[[i]], pars = param, probs = c(0.025, 0.5, 0.975), digits_summary = 4)$summary

                 mean      se_mean          sd       2.5%        50%
intercept  0.02279584 0.0005675383 0.004899672 0.01377460 0.02295105
lscale[1]  0.16321279 0.0012264089 0.016537237 0.12621907 0.16452300
lscale[2]  0.22884513 0.0024885582 0.021755035 0.18737257 0.22838250
lscale[3]  1.78115501 0.0484845953 0.621830253 0.73027560 1.74175500
gpscale[1] 0.62287706 0.0104260516 0.142880281 0.43168012 0.59825100
gpscale[2] 0.29488222 0.0048951827 0.055301127 0.21142265 0.28873300
gpscale[3] 3.63075626 0.2139917502 2.737805014 0.72485948 2.59414000
noise      0.29303347 0.0001938712 0.002414114 0.28847757 0.29308650
sigma_f4   0.01280845 0.0003352905 0.002444182 0.00896321 0.01251480
                 97.5%     n_eff      Rhat
intercept   0.03294064  74.53218 1.0714969
lscale[1]   0.19295540 181.82570 1.0068069
lscale[2]   0.27628600  76.42298 1.0651080
lscale[3]   3.06536775 164.48873 1.0029784
gpscale[1]  0.95707893 187.80400 0.9938507
gpscale[2]  0.43470313 127.62335 1.0243951
gpscale[3] 10.83963750 163.68579 1.0195254
noise       0.29736762 155.05609 1.0005674
sigma_f4    0.01857094  53.14037 1.0729330

Simulation chains for the variables after warmup

traceplot(fit[[i]], pars = param, include = TRUE, unconstrain = FALSE, inc_warmup = FALSE, window = NULL, nrow = NULL, ncol = NULL)

Plot of the mean posterior function \(f_1\)

Plot of the mean posterior function \(f_2\)

Plot of the mean posterior function \(f_1 + f_2\)

Plot of the residuals

Model evaluation

Residuals

f <- summary(fit[[i]], pars = c("f"))$summary[,1]
res[[i]] <- standata[[i]]$y - f

ggplot(as.data.frame(res[[i]]), aes(res[[i]]))  +
  geom_histogram(color = 'white') +
  theme_classic()

Root mean squared error

rmse[i] <- sqrt(mean(res[[i]]^2))
rmse[i]

[1] 0.2895321

Bayesian \(R^2\) (Coefficient of determination)

noise <- as.matrix(fit[[i]], pars = c("noise"))
sd_f <- apply(as.matrix(fit[[i]], pars = c("f")), 1, sd)
R2 <- sd_f^2/(sd_f^2 + noise^2)

ggplot(as.data.frame(R2), aes(R2))  +
  geom_histogram(color = 'white') +
  theme_classic()

Median of \(R^2\)

eR2[i] <- median(R2)
eR2[i]

[1] 0.9140058

Log predictive density (lpd)

lpd <- summary(fit[[i]], pars = c("log_lik"))$summary[,1]

ggplot(as.data.frame(lpd), aes(lpd))  +
  geom_histogram(color = 'white') +
  theme_classic()

Median of lpd

elpd[i] <- median(lpd)
elpd[i]

[1] 0.1497226

Diagnosis

Smooth trend function \(f_1\)

1. Estimated lengthscale \(\hat{l}\)

l_f1_hat[i] <- round(summary(fit[[i]], pars = "lscale[1]")$summary[,1], 2)
l_f1_hat[i]

[1] 0.16

2. Check whether \(\hat{l}\) is equal to or greater than l_f1 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f1_hat[i] >= l_f1[i]

[1] TRUE

Yearly periodic effects function \(f_2\)

1. Estimated lengthscale \(\hat{l}\)

l_f2_hat[i] <- round(summary(fit[[i]], pars = "lscale[2]")$summary[,1], 2)
l_f2_hat[i]

[1] 0.23

2. Check whether \(\hat{l}\) is equal to or greater than l_f2 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f2_hat[i] >= l_f2[i]

[1] FALSE

Weekly periodic effects function \(f_3\)

1. Estimated lengthscale \(\hat{l}\)

l_f3_hat[i] <- round(summary(fit[[i]], pars = "lscale[3]")$summary[,1], 2)
l_f3_hat[i]

[1] 1.78

2. Check whether \(\hat{l}\) is equal to or greater than l_f3 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f3_hat[i] >= l_f3[i]

[1] TRUE

Summary table

iter	GP_func	l	c	m	l_hat	l_hat > l	rmse	R2	elpd
5	f1	0.16	1.5	32	0.16	TRUE	0.29	0.914	0.15
5	f2	0.25	NA	15	0.23	FALSE	0.29	0.914	0.15
5	f3	1.64	NA	5	1.78	TRUE	0.29	0.914	0.15

Iteration 6

Iteration index

i <- 6

Setting \(m\), \(l\) and \(c\)

Smooth trend function \(f_1\)

1. Updating \(l\) with its estimate \(\hat{l}\) in the fifth iteration

l_f1[i] <- l_f1_hat[i-1]
l_f1[i]

[1] 0.16

2. The boundary factor \(c\) has to fulfill \(c \, \geq \, 0.97 + 1.45 \, \frac{l}{S}\)

0.97 + 1.45*l_f1[i]/S

[1] 1.103973

Then, \(c\) is set to

c_f1[i] <- 1.5
c_f1[i]

[1] 1.5

3. The number of basis functions \(m\) as a function of \(l\): \(\; m= 1.73 \, \frac{c}{(l/S)^{1.05}}\)

if(1.73 * c_f1[i]/(l_f1[i]/S)^1.05 < 5) m_f1[i] <- 5 else m_f1[i] <- ceiling(1.73 * c_f1[i]/(l_f1[i]/S)^1.05)
m_f1[i]

[1] 32

Yearly periodic effects function \(f_2\)

1. Updating \(l\) with its estimate \(\hat{l}\) in the fifth iteration

l_f2[i] <- l_f2_hat[i-1]
l_f2[i]

[1] 0.23

2. Number of basis functions \(m\) as a function of \(l\): \(\; m= 3.31 \, \frac{1}{l^{1.05}}\)

if(3.31 * 1/l_f2[i]^1.05 < 5) m_f2[i] <- 5 else m_f2[i] <- ceiling(3.31 * 1/l_f2[i]^1.05)
m_f2[i]

[1] 16

Weekly periodic effects function \(f_3\)

1. Updating \(l\) with its estimate \(\hat{l}\) in the fifth iteration

l_f3[i] <- l_f3_hat[i-1]
l_f3[i]

[1] 1.78

2. Number of basis functions \(m\) as a function of \(l\): \(\; m= 3.31 \, \frac{1}{l^{1.05}}\)

if(3.31 * 1/l_f3[i]^1.05 < 5) m_f3[i] <- 5 else m_f3[i] <- ceiling(3.31 * 1/l_f3[i]^1.05)
m_f3[i]

[1] 5

Data to Stan

standata[[i]] <- list(M_f1= m_f1[i], 
                 c_f1= c_f1[i], 
                 J_f2= m_f2[i], 
                 J_f3= m_f3[i], 
                 x= x[,1], 
                 y= y[,1], 
                 N= length(x), 
                 period_year= period_year, 
                 period_week= period_week,
                 day_of_year2= day_of_year2
)
str(standata[[i]])

List of 10
 $ M_f1        : num 32
 $ c_f1        : num 1.5
 $ J_f2        : num 16
 $ J_f3        : num 5
 $ x           : num [1:7305] -1.73 -1.73 -1.73 -1.73 -1.73 ...
 $ y           : num [1:7305] -1.0316 -0.5739 -0.0949 -0.6111 -1.1168 ...
 $ N           : int 7305
 $ period_year : num 0.173
 $ period_week : num 0.00332
 $ day_of_year2: int [1:7305] 1 2 3 4 5 6 7 8 9 10 ...

Model fitting

Compiling the model

# birthday_mod <- cmdstanr::cmdstan_model(stan_file = "stancode_def.stan")

Model sampling

# load("fit.rData")

# fit[[i]] <- birthday_mod$sample(data= standata[[i]], iter_warmup=200, iter_sampling=200, chains=4, thin=4, init=0.5, adapt_delta=0.9, save_warmup=FALSE)
# fit[[i]] <- rstan::read_stan_csv(fit[[i]]$output_files())

# save(fit, file="fit.rData")

Summaries of variable estimates

param = c("intercept","lscale[1]","lscale[2]","lscale[3]","gpscale[1]","gpscale[2]","gpscale[3]","noise","sigma_f4")
summary(fit[[i]], pars = param, probs = c(0.025, 0.5, 0.975), digits_summary = 4)$summary

                 mean      se_mean          sd        2.5%        50%
intercept  0.02294053 0.0004648112 0.005550149 0.013452927 0.02243810
lscale[1]  0.16117096 0.0012463943 0.016428119 0.124791725 0.16174750
lscale[2]  0.23273374 0.0017141776 0.020520817 0.196417325 0.23225950
lscale[3]  1.64564855 0.0541690077 0.644079224 0.618704450 1.61698500
gpscale[1] 0.58483434 0.0120735964 0.131891179 0.411753125 0.56481850
gpscale[2] 0.30437833 0.0042693583 0.051882474 0.223969375 0.29740400
gpscale[3] 3.16641534 0.2283539921 2.821447927 0.656188800 2.26845000
noise      0.29283316 0.0001855130 0.002660376 0.287663800 0.29305100
sigma_f4   0.01399103 0.0004748107 0.002561988 0.009126358 0.01402015
                97.5%     n_eff      Rhat
intercept  0.03383408 142.57918 0.9886644
lscale[1]  0.19031025 173.72597 1.0078613
lscale[2]  0.27478578 143.31039 1.0149307
lscale[3]  3.11322025 141.37637 1.0113828
gpscale[1] 0.88167815 119.33235 1.0145168
gpscale[2] 0.43783017 147.67812 1.0041472
gpscale[3] 9.64981850 152.66056 1.0196156
noise      0.29808180 205.65415 0.9916026
sigma_f4   0.01883736  29.11475 1.1056103

Simulation chains for the variables after warmup

traceplot(fit[[i]], pars = param, include = TRUE, unconstrain = FALSE, inc_warmup = FALSE, window = NULL, nrow = NULL, ncol = NULL)

Plot of the mean posterior function \(f_1\)

Plot of the mean posterior function \(f_2\)

Plot of the mean posterior function \(f_1 + f_2\)

Plot of the residuals

Model evaluation

Residuals

f <- summary(fit[[i]], pars = c("f"))$summary[,1]
res[[i]] <- standata[[i]]$y - f

ggplot(as.data.frame(res[[i]]), aes(res[[i]]))  +
  geom_histogram(color = 'white') +
  theme_classic()

Root mean squared error

rmse[i] <- sqrt(mean(res[[i]]^2))
rmse[i]

[1] 0.2893593

Bayesian \(R^2\) (Coefficient of determination)

noise <- as.matrix(fit[[i]], pars = c("noise"))
sd_f <- apply(as.matrix(fit[[i]], pars = c("f")), 1, sd)
R2 <- sd_f^2/(sd_f^2 + noise^2)

ggplot(as.data.frame(R2), aes(R2))  +
  geom_histogram(color = 'white') +
  theme_classic()

Median of \(R^2\)

eR2[i] <- median(R2)
eR2[i]

[1] 0.9142495

Log predictive density (lpd)

lpd <- summary(fit[[i]], pars = c("log_lik"))$summary[,1]

ggplot(as.data.frame(lpd), aes(lpd))  +
  geom_histogram(color = 'white') +
  theme_classic()

Median of lpd

elpd[i] <- median(lpd)
elpd[i]

[1] 0.1472021

Diagnosis

Smooth trend function \(f_1\)

1. Estimated lengthscale \(\hat{l}\)

l_f1_hat[i] <- round(summary(fit[[i]], pars = "lscale[1]")$summary[,1], 2)
l_f1_hat[i]

[1] 0.16

2. Check whether \(\hat{l}\) is equal to or greater than l_f1 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f1_hat[i] >= l_f1[i]

[1] TRUE

Yearly periodic effects function \(f_2\)

1. Estimated lengthscale \(\hat{l}\)

l_f2_hat[i] <- round(summary(fit[[i]], pars = "lscale[2]")$summary[,1], 2)
l_f2_hat[i]

[1] 0.23

2. Check whether \(\hat{l}\) is equal to or greater than l_f2 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f2_hat[i] >= l_f2[i]

[1] TRUE

Weekly periodic effects function \(f_3\)

1. Estimated lengthscale \(\hat{l}\)

l_f3_hat[i] <- round(summary(fit[[i]], pars = "lscale[3]")$summary[,1], 2)
l_f3_hat[i]

[1] 1.65

2. Check whether \(\hat{l}\) is equal to or greater than l_f3 (the minimum \(l\) that can be accurately fitted determined by the \(m\) used)

l_f3_hat[i] >= l_f3[i]

[1] FALSE

Summary table

iter	GP_func	l	c	m	l_hat	l_hat > l	rmse	R2	elpd
6	f1	0.16	1.5	32	0.16	TRUE	0.289	0.914	0.147
6	f2	0.23	NA	16	0.23	TRUE	0.289	0.914	0.147
6	f3	1.78	NA	5	1.65	FALSE	0.289	0.914	0.147

Global summary of the diagnosis (with all iterations)

	iter	GP_func	l	c	m	l_hat	l_hat > l	rmse	R2	elpd
1	1	f1	0.52	1.5	10	0.24	FALSE	0.306	0.904	0.079
4	2	f1	0.24	1.5	21	0.20	FALSE	0.295	0.910	0.118
7	3	f1	0.20	1.5	26	0.13	FALSE	0.293	0.912	0.132
10	4	f1	0.13	1.5	40	0.16	TRUE	0.289	0.914	0.149
13	5	f1	0.16	1.5	32	0.16	TRUE	0.290	0.914	0.150
16	6	f1	0.16	1.5	32	0.16	TRUE	0.289	0.914	0.147

	iter	GP_func	l	c	m	l_hat	l_hat > l	rmse	R2	elpd
2	1	f2	0.50	NA	7	0.44	FALSE	0.306	0.904	0.079
5	2	f2	0.44	NA	8	0.35	FALSE	0.295	0.910	0.118
8	3	f2	0.35	NA	10	0.29	FALSE	0.293	0.912	0.132
11	4	f2	0.29	NA	13	0.25	FALSE	0.289	0.914	0.149
14	5	f2	0.25	NA	15	0.23	FALSE	0.290	0.914	0.150
17	6	f2	0.23	NA	16	0.23	TRUE	0.289	0.914	0.147

	iter	GP_func	l	c	m	l_hat	l_hat > l	rmse	R2	elpd
3	1	f3	0.50	NA	7	1.59	TRUE	0.306	0.904	0.079
6	2	f3	1.59	NA	5	1.70	TRUE	0.295	0.910	0.118
9	3	f3	1.70	NA	5	1.57	FALSE	0.293	0.912	0.132
12	4	f3	1.57	NA	5	1.64	TRUE	0.289	0.914	0.149
15	5	f3	1.64	NA	5	1.78	TRUE	0.290	0.914	0.150
18	6	f3	1.78	NA	5	1.65	FALSE	0.289	0.914	0.147

Plots of the mean posterior functions from the last iteration

Extracting estimated function components \(f_1\), \(f_2\), \(f_3\) and \(f_4\)

f <- summary(fit[[i]], pars = c("f"), probs = c(0.025, 0.5, 0.975))$summary
intercept <- summary(fit[[i]], pars = c("intercept"), probs = c(0.025, 0.5, 0.975))$summary
f1 <- summary(fit[[i]], pars = c("f1"), probs = c(0.025, 0.5, 0.975))$summary
f2 <- summary(fit[[i]], pars = c("f2"), probs = c(0.025, 0.5, 0.975))$summary
f3 <- summary(fit[[i]], pars = c("f3"), probs = c(0.025, 0.5, 0.975))$summary
f4 <- summary(fit[[i]], pars = c("f4"), probs = c(0.025, 0.5, 0.975))$summary

Plot of only one year, the year 1972

Plot of all the years

Plot of only the first month

Plot of the first four years